home *** CD-ROM | disk | FTP | other *** search
/ Magnum One / Magnum One (Mid-American Digital) (Disc Manufacturing).iso / d12 / snip9_91.arc / RAND1.C < prev    next >
C/C++ Source or Header  |  1991-09-17  |  8KB  |  237 lines
  1. /************************************************************************
  2.  This random number generator originally appeared in "Toward a Universal
  3.  Random Number Generator" by George Marsaglia and Arif Zaman.
  4.  Florida State University Report: FSU-SCRI-87-50 (1987)
  5.  
  6.  It was later modified by F. James and published in "A Review of Pseudo-
  7.  random Number Generators"
  8.  
  9.  Converted from FORTRAN to C by Phil Linttell, James F. Hickling
  10.  Management Consultants Ltd, Aug. 14, 1989.
  11.  
  12.  THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
  13.        (However, a newly discovered technique can yield
  14.          a period of 10^600. But that is still in the development stage.)
  15.  
  16.  It passes ALL of the tests for random number generators and has a period
  17.    of 2^144, is completely portable (gives bit identical results on all
  18.    machines with at least 24-bit mantissas in the floating point
  19.    representation).
  20.  
  21.  The algorithm is a combination of a Fibonacci sequence (with lags of 97
  22.    and 33, and operation "subtraction plus one, modulo one") and an
  23.    "arithmetic sequence" (using subtraction).
  24.  
  25.  On a Vax 11/780, this random number generator can produce a number in
  26.     13 microseconds.
  27. ************************************************************************/
  28.  
  29. #include <stdio.h>
  30. #include <stdlib.h>
  31. #include <math.h>
  32.  
  33. #define TRUE    1
  34. #define FALSE   0
  35.  
  36. float u[97], c, cd, cm;
  37. int i97, j97, test;
  38.  
  39. int rmarin(int ij, int kl);
  40. int ranmar(float rvec[], int len);
  41.  
  42.  
  43. int main()
  44. {
  45.  
  46.         float temp[100];
  47.         int i;
  48.         int ij, kl, len;
  49.  
  50.         /*These are the seeds needed to produce the test case results*/
  51.  
  52.         ij = 1802;
  53.         kl = 9373;
  54.  
  55.         /*Do the initialization*/
  56.  
  57.         if (1 == rmarin(ij,kl))
  58.                 return 1;
  59.  
  60.         /*Generate 20000 random numbers*/
  61.  
  62.         len = 100;
  63.         for ( i=0; i<=199 ; i++)
  64.                 if (1 == ranmar(temp, len))
  65.                         return 1;
  66.     
  67.         /*If the random number generator is working properly,
  68.           the next six random numbers should be:
  69.  
  70.             6533892.0  14220222.0   7275067.0
  71.             6172232.0   8354498.0  10633180.0
  72.         */
  73.  
  74.         len = 6;
  75.         if (1 == ranmar(temp, len))
  76.                 return 1;
  77.  
  78.         for ( i=0; i<=5; i++)
  79.                 printf("%12.1f\n",4096.0*4096.0*temp[i]);
  80.  
  81.         return NULL;
  82. }
  83.  
  84.  
  85. /************************************************************************
  86.  This is the initialization routine for the random number generator RANMAR()
  87.  NOTE: The seed variables can have values between:    0 <= IJ <= 31328
  88.                                                       0 <= KL <= 30081
  89.  The random number sequences created by these two seeds are of sufficient
  90.  length to complete an entire calculation with. For example, if several
  91.  different groups are working on different parts of the same calculation,
  92.  each group could be assigned its own IJ seed. This would leave each group
  93.  with 30000 choices for the second seed. That is to say, this random
  94.  number generator can create 900 million different subsequences -- with
  95.  each subsequence having a length of approximately 10^30.
  96.  
  97.  Use IJ = 1802 & KL = 9373 to test the random number generator. The
  98.  subroutine RANMAR should be used to generate 20000 random numbers.
  99.  Then display the next six random numbers generated multiplied by 4096*4096
  100.  If the random number generator is working properly, the random numbers
  101.  should be:
  102.            6533892.0  14220222.0   7275067.0
  103.            6172232.0   8354498.0  10633180.0
  104. ************************************************************************/
  105.  
  106. int rmarin(int ij, int kl)
  107. {
  108.  
  109.         float s, t;
  110.         int i, j, k, l, m;
  111.         int ii, jj;
  112.  
  113.         /* Change FALSE to TRUE in the next statement to test the
  114.            random routine.*/
  115.  
  116.         test = TRUE;
  117.  
  118.         if ( ( ij < 0 || ij > 31328 ) ||
  119.                 ( kl < 0 || kl > 30081 ) )
  120.         {
  121.                 printf ("RMARIN: The first random number seed must have a "
  122.                         "value between 0 and 31328\n");
  123.                 printf ("        The second random number seed must have a "
  124.                         "value between 0 and 30081");
  125.                 return 1;
  126.         }
  127.  
  128.         i = fmod(ij/177.0, 177.0) + 2;
  129.         j = fmod(ij      , 177.0) + 2;
  130.         k = fmod(kl/169.0, 178.0) + 1;
  131.         l = fmod(kl      , 169.0);
  132.  
  133.         for ( ii=0; ii<=96; ii++ )
  134.         {
  135.                 s = 0.0;
  136.                 t = 0.5;
  137.                 for ( jj=0; jj<=23; jj++ )
  138.                 {
  139.                         m = fmod( fmod(i*j,179.0)*k , 179.0 );
  140.                         i = j;
  141.                         j = k;
  142.                         k = m;
  143.                         l = fmod( 53.0*l+1.0 , 169.0 );
  144.                         if ( fmod(l*m,64.0) >= 32)
  145.                                 s = s + t;
  146.                         t = 0.5 * t;
  147.                 }
  148.                 u[ii] = s;
  149.         }
  150.  
  151.         c  =   362436.0 / 16777216.0;
  152.         cd =  7654321.0 / 16777216.0;
  153.         cm = 16777213.0 / 16777216.0;
  154.  
  155.         i97 = 96;
  156.         j97 = 32;
  157.  
  158.         test = TRUE;
  159.  
  160.         return NULL;
  161. }
  162.  
  163. int ranmar(float rvec[], int len)
  164. {
  165.         float uni;
  166.         int ivec;
  167.  
  168.         if ( !test )
  169.         {
  170.                 printf ("RANMAR: Call the initialization routine (RMARIN) "
  171.                         "before calling RANMAR.\n");
  172.                 return 1;
  173.         }
  174.  
  175.         for ( ivec=0; ivec < len; ivec++)
  176.         {
  177.                 uni = u[i97] - u[j97];
  178.                 if ( uni < 0.0 )
  179.                         uni = uni + 1.0;
  180.                 u[i97] = uni;
  181.                 i97--;
  182.                 if ( i97 < 0 )
  183.                         i97 = 96;
  184.                 j97--;
  185.                 if ( j97 < 0 )
  186.                         j97 = 96;
  187.                 c = c - cd;
  188.                 if ( c < 0.0 )
  189.                         c = c + cm;
  190.                 uni = uni - c;
  191.                 if ( uni < 0.0 )
  192.                         uni = uni + 1.0;
  193.                 rvec[ivec] = uni;
  194.         }
  195.         return NULL;
  196. }
  197.  
  198. /* I use the following procedure in TC to generate seeds:
  199.  
  200.   The sow() procedure calculates two seeds for use with the random number
  201.   generator from the system clock.  I decided how to do this myself, and
  202.   I am sure that there must be better ways to select seeds; hopefully,
  203.   however, this is good enough.  The first seed is calculated from the values
  204.   for second, minute, hour, and year-day; weighted with the second most
  205.   significant and year-day least significant.  The second seed weights the
  206.   values in reverse.
  207. */
  208.  
  209. void sow( seed1, seed2 )
  210. int *seed1, *seed2;
  211. {
  212.         struct tm *tm_now;
  213.         float s_sig, s_insig, maxs_sig, maxs_insig;
  214.         long secs_now;
  215.         int s, m, h, d, s1, s2;
  216.  
  217.         time(&secs_now);
  218.         tm_now = localtime(&secs_now);
  219.  
  220.         s = tm_now->tm_sec + 1;
  221.         m = tm_now->tm_min + 1;
  222.         h = tm_now->tm_hour + 1;
  223.         d = tm_now->tm_yday + 1;
  224.  
  225.         maxs_sig   = 60.0 + 60.0/60.0 + 24.0/60.0/60.0 + 366.0/24.0/60.0/60.0;
  226.         maxs_insig = 60.0 + 60.0*60.0 + 24.0*60.0*60.0 + 366.0*24.0*60.0*60.0;
  227.  
  228.         s_sig      = s + m/60.0 + h/60.0/60.0 + d/24.0/60.0/60.0;
  229.         s_insig    = s + m*60.0 + h*60.0*60.0 + d*24.0*60.0*60.0;
  230.  
  231.         s1 = s_sig   / maxs_sig   * 31328.0;
  232.         s2 = s_insig / maxs_insig * 30081.0;
  233.  
  234.         *seed1 = s1;
  235.         *seed2 = s2;
  236. }
  237.